% Read the U and V data from Neetu's Chile simulation netCDF file
% Interpolate the data to 0.125 intervals
% Fourier transform and save the data 

%model27feb10.nc

%  dimensions/  Group  
%  LONGITUDEN560_N160 3.32 K Dataset 401 values -159.9170  -59.9170  at 0.25 sep
%  LATITUDE81_321 2.07 K Dataset  241 values  -59.9170 to 0.0830 at 0.025 deg
%  TIME 5.81 K Dataset  
%  SL 265.43 M Dataset  
%  UUT 265.43 M Dataset  
%  VVT 265.43 M Dataset  
%  history 180 B Attribute 

nc_fname = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\model27feb10.nc';
ncid = netcdf.open(nc_fname,'NC_WRITE');
sea_level_id = netcdf.inqVarID(ncid,'SL');
long_id = netcdf.inqVarID(ncid,'LONGITUDEN560_N160');
lat_id = netcdf.inqVarID(ncid,'LATITUDE81_321');
time_id = netcdf.inqVarID(ncid,'TIME');
U_id = netcdf.inqVarID(ncid,'UUT');
V_id = netcdf.inqVarID(ncid,'VVT');
long_data = netcdf.getVar(ncid,long_id,0,401);
lat_data = netcdf.getVar(ncid,lat_id,0,241);
time_data = netcdf.getVar(ncid,time_id,0,720);

% sea_level_data =  netcdf.getVar(ncid,3,[0 0 0],[401 241 720]);
% fillValue = min(min(min(sea_level_data))); 
% L = sea_level_data == fillValue;
% sea_level_data(L) = NaN;

[OX,OY] = meshgrid([-158.875:0.25:-60.125]',[-58.875:0.25:0.125]'); %why these values ? = see the limit of lat1 and lon1


U_data =  netcdf.getVar(ncid,U_id,[0 0 0],[401 241 720]);
fillValue = max(max(max(U_data)));
L = U_data == fillValue;
U_data(L) = NaN;


for i = 1:360,
    IU_n_1(:,:,i) = single(interp2(long_data,lat_data,squeeze(U_data(:,:,i))',OX,OY,'nearest'));
end;

save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IU_n1.mat IU_n_1;
clear IU_n_1;

for i = 361:720,
    IU_n_2(:,:,i) = single(interp2(long_data,lat_data,squeeze(U_data(:,:,i))',OX,OY,'nearest'));
end;

save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IU_n2.mat IU_n_2;
clear IU_n_2;
clear U_data;

V_data =  netcdf.getVar(ncid,V_id,[0 0 0],[401 241 720]);
fillValue = max(max(max(V_data)));
L = V_data == fillValue;
V_data(L) = NaN;


for i = 1:360,
    IV_n_1(:,:,i) = single(interp2(long_data,lat_data,squeeze(V_data(:,:,i))',OX,OY,'nearest'));
end;

save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IV_n1.mat IV_n_1;
clear IV_n_1;

for i = 361:720,
    IV_n_2(:,:,i) = single(interp2(long_data,lat_data,squeeze(V_data(:,:,i))',OX,OY,'nearest'));
end;

save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IV_n2.mat IV_n_2;
clear IV_n_2;


%% fft
load d:\manoj\tsunami\2010_Chile_Easter_Island_data\modeling\IV;
IV_spec = fft(IV_n_2,[],3);
save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IV_spec IV_spec;
load d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IU;
IU_spec = fft(IU_n_2,[],3);
save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\IU_spec IU_spec;

prd = 1./(fftfrq(720,1/(60/(24*3600))));
prd = [4000. prd]; % The mean flow period -> 4000 days (it should have been DC but 4000 days OK)!


%% Bz
bz_025deg = load('c:\manoj\temp\ts_simulation.txt');
bz = reshape(bz_025deg(:,5),[1440,720]);
save d:\manoj\tsunami\2010_Chile_Easter_Island_data\Velocity\bz_wmm2010_025deg bz;

%% take the 0.25 conductance from the file, remove the unwanted latitudes and put back the file (model)
conductance = load('d:\manoj\ocean\x3dg\conductance025.txt');
con1 = reshape(conductance',[1440,720]);
%conduct = con1(:,360:596);%remove the unwanted latitudes
conduct = con1(:,420:519);%remove the unwanted latitudes (420 corresponding to -14.875 lat)
%con1 = reshape(conduct,[18,18960])';% put back into the same array as in the file
con1 = reshape(conduct,[18,8000])';% put back into the same array as in the file

fid = fopen('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\latitude_limited_conductance.txt','wt');
    for j = 1:8000,
        fprintf(fid,'%9.2f',con1(j,:));
        fprintf(fid,'\n');
    end;

    %% Print latitude array
    
for i = 1:100,
fprintf('%5.2f',0.25);
if mod(i,25)  == 0,
fprintf('\n');
end;
end;
    
%%
data_ur = zeros([237,1440]);
data_ui = zeros([237,1440]);
data_vr = zeros([237,1440]);
data_vi = zeros([237,1440]);
prd = 1./(fftfrq(237,1/(60/(24*3600))));
prd = [4000. prd]; % The mean flow period -> 4000 days (it should have been DC but 4000 days OK)!
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\bz_wmm2010_025deg bz;
bz = rot90(bz);
Bz = bz(360:596,:);
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\IU_spec;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\IV_spec;

i = 2;
%data_ur(360:596,805:1200) = flipud(real(squeeze(IU_spec(:,:,i))));
data_ur(:,805:1200) = flipud(real(squeeze(IU_spec(:,:,i))));
data_ui(:,805:1200) = flipud(imag(squeeze(IU_spec(:,:,i))));
data_vr(:,805:1200) = flipud(real(squeeze(IV_spec(:,:,i))));
data_vi(:,805:1200) = flipud(imag(squeeze(IV_spec(:,:,i))));
    
%Note: 26July06. The original grid Neetu (new comp) gave was
%
    
CUr = -data_vr.*Bz.*1e-9.*3.2;
CUi = -data_vi.*Bz.*1e-9.*3.2;
CVr = -data_ur.*Bz.*1e-9.*3.2;
CVi = -data_ui.*Bz.*1e-9.*3.2;
%%
%fprintf(fp,'%s\n',['    ' sprintf('%10.5f',prd(i)) '     ' sprintf('tsu%d',i) ' 1x1_mantle_k3.model ' sprintf('tsu%d',i) '.source' ' 1']);
fid = fopen('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\test.source','wt');
fid1 = fopen('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\header_deg0.25.txt','rt');
%%
%%%the following format for alexei's x3dg is tested ok 20 oct 05
           for ii = 1:34,
                s=fgetl(fid1);
                fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:237,
    for jn = 1:1440,
       fprintf(fid,'%12.4e\n',CUr(in,jn));
   end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:237,
    for jn = 1:1440,
        fprintf(fid,'%12.4e\n',CUi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
 for in = 1:237,
    for jn = 1:1440,       
        fprintf(fid,'%12.4e\n',CVr(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
  
 for in = 1:237,
    for jn = 1:1440,       
        fprintf(fid,'%12.4e\n',CVi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:37,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
            fclose(fid);
            fclose(fid1);

fclose all;

%%

data_ur = zeros([100,1440]);
data_ui = zeros([100,1440]);
data_vr = zeros([100,1440]);
data_vi = zeros([100,1440]);
prd = 1./(fftfrq(720,1/(60/(24*3600))));
prd = [4000. prd]; % The mean flow period -> 4000 days (it should have been DC but 4000 days OK)!
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\bz_wmm2010_025deg bz;
bz = rot90(bz);
Bz = bz(420:519,:);
%find((89.875:-0.25:-89.875)==0.125) = 360 (cell number)
%cell number for north lat limit -14.875 is 420, same for lat -39.625 is
%519
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\IU_spec;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\IV_spec;
%%
fid_aproject = fopen('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\tsunami_a.project','at');
%301
for i = 301 : length(prd),


%data_ur(i,j) i for lat and j for long lat limits are set by the 
% definitions above, long indexes correspond to the long value limits of
% the IU and IV data.  Long limits are [-158.875:0.25:-60.125]',[-58.875:0.25:0.125]'
% find(0.125:0.25:359.875==(360-158.875)) = 805
% the lat limits of IU and IV are [-58.875:0.25:0.125]'
% so to take only date between -14.875 AND -39.625
% No ifort please, use g77 to compile x3dg_dr_reg.f
data_ur(:,805:1200) = flipud(real(squeeze(IU_spec(78:177,:,i))));
data_ui(:,805:1200) = flipud(imag(squeeze(IU_spec(78:177,:,i))));
data_vr(:,805:1200) = flipud(real(squeeze(IV_spec(78:177,:,i))));
data_vi(:,805:1200) = flipud(imag(squeeze(IV_spec(78:177,:,i))));
    
%Note: 26July06. The original grid Neetu (new comp) gave was
%
    
CUr = -data_vr.*Bz.*1e-9.*3.2;
CUi = -data_vi.*Bz.*1e-9.*3.2;
CVr = -data_ur.*Bz.*1e-9.*3.2;
CVi = -data_ui.*Bz.*1e-9.*3.2;

array_dim = size(CUr);

fprintf(fid_aproject,'%s\n',['    ' sprintf('%10.5f',prd(i)) '     ' sprintf('tsu%d',i) ' tsunami25.model ' sprintf('tsu%d',i) '.source' ' 1']);
display(['    ' sprintf('%10.5f',prd(i)) '     ' sprintf('tsu%d',i) ' tsunami25.model ' sprintf('tsu%d',i) '.source' ' 1']);
fid = fopen(['D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\sourcefiles\' sprintf('tsu%d',i) '.source'],'wt');
fid1 = fopen('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\header_deg0.25.txt','rt');

%%%the following format for alexei's x3dg is tested ok 20 oct 05
           for ii = 1:34,
                s=fgetl(fid1);
                fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:array_dim(1),
    for jn = 1:array_dim(2),
       fprintf(fid,'%12.4e\n',CUr(in,jn));
   end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:array_dim(1),
    for jn = 1:array_dim(2),
        fprintf(fid,'%12.4e\n',CUi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
 for in = 1:array_dim(1),
    for jn = 1:array_dim(2),       
        fprintf(fid,'%12.4e\n',CVr(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
  
 for in = 1:array_dim(1),
    for jn = 1:array_dim(2),       
        fprintf(fid,'%12.4e\n',CVi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:37,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
            fclose(fid);
            fclose(fid1);



end

fclose all;
